Linear Mixed Models
Single Random Effect
data(reeds, package="brinla")
summary(reeds)
site nitrogen
A:5 Min. :2.35
B:5 1st Qu.:2.62
C:5 Median :3.06
Mean :3.03
3rd Qu.:3.27
Max. :3.93
library(lme4)
mmod <- lmer(nitrogen ~ 1+(1|site), reeds)
summary(mmod)
Linear mixed model fit by REML ['lmerMod']
Formula: nitrogen ~ 1 + (1 | site)
Data: reeds
REML criterion at convergence: 13
Scaled residuals:
Min 1Q Median 3Q Max
-1.221 -0.754 -0.263 0.914 1.612
Random effects:
Groups Name Variance Std.Dev.
site (Intercept) 0.1872 0.433
Residual 0.0855 0.292
Number of obs: 15, groups: site, 3
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.029 0.261 11.6
library(INLA)
formula <- nitrogen ~ 1 + f(site, model="iid")
imod <- inla(formula,family="gaussian", data = reeds)
summary(imod)
Call:
"inla(formula = formula, family = \"gaussian\", data = reeds)"
Time used:
Pre-processing Running inla Post-processing Total
0.9781 0.1516 0.2377 1.3674
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 3.0293 0.1402 2.7587 3.0293 3.3004 3.0293 1e-04
Random effects:
Name Model
site IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 13.25 5.203 5.471 12.49 25.59 10.940
Precision for site 20.35 25.120 1.933 12.81 84.76 5.093
Expected number of effective parameters(std dev): 1.243(0.5604)
Number of equivalent replicates : 12.07
Marginal log-Likelihood: -20.88
invsqrt <- function(x) 1/sqrt(x)
sdt <- invsqrt(imod$summary.hyperpar[,-2])
row.names(sdt) <- c("SD of epsilson","SD of site")
sdt
mean 0.025quant 0.5quant 0.975quant mode
SD of epsilson 0.27471 0.42753 0.28294 0.19767 0.30233
SD of site 0.22168 0.71916 0.27942 0.10862 0.44309
prec.site <- imod$marginals.hyperpar$"Precision for site"
prec.epsilon <- imod$marginals.hyperpar$"Precision for the Gaussian observations"
c(epsilon=inla.emarginal(invsqrt,prec.epsilon),site=inla.emarginal(invsqrt,prec.site))
epsilon site
0.29086 0.31354
sigma.site <- inla.tmarginal(invsqrt, prec.site)
sigma.epsilon <- inla.tmarginal(invsqrt, prec.epsilon)
c(epsilon=inla.mmarginal(sigma.epsilon),site=inla.mmarginal(sigma.site))
epsilon site
0.26648 0.22175
inla.contrib.sd(imod)$hyper
mean sd 2.5% 97.5%
sd for the Gaussian observations 0.28960 0.058561 0.19695 0.42688
sd for site 0.31705 0.165926 0.11014 0.76824
sampvars <- 1/inla.hyperpar.sample(1000,imod)
sampicc <- sampvars[,2]/(rowSums(sampvars))
quantile(sampicc, c(0.025,0.5,0.975))
2.5% 50% 97.5%
0.086604 0.495696 0.898860
library(brinla)
bri.hyperpar.summary(imod)
mean sd q0.025 q0.5 q0.975 mode
SD for the Gaussian observations 0.29079 0.05833 0.19817 0.28283 0.42608 0.26647
SD for site 0.31313 0.15740 0.10893 0.27909 0.71371 0.22176
alpha <- data.frame(imod$marginals.fixed[[1]])
library(ggplot2)
ggplot(alpha, aes(x,y)) + geom_line() + xlim(2,4) + geom_vline(xintercept = c(2.66, 3.40))+xlab("nitrogen")+ylab("density")
x <- seq(0,1,len=100)
d1 <- inla.dmarginal(x,sigma.site)
d2 <- inla.dmarginal(x,sigma.epsilon)
rdf <- data.frame(nitrogen=c(x,x),sigma=gl(2,100,labels=c("site","epsilon")),density=c(d1,d2))
ggplot(rdf, aes(x=nitrogen, y=density, linetype=sigma))+geom_line()
bri.hyperpar.plot(imod)
sdres <- sd(reeds$nitrogen)
pcprior <- list(prec = list(prior="pc.prec", param = c(3*sdres,0.01)))
formula <- nitrogen ~ f(site, model="iid", hyper = pcprior)
pmod <- inla(formula, family="gaussian", data=reeds)
pmod$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 3.0293 0.30703 2.4026 3.0293 3.6569 3.0293 2.2118e-07
bri.hyperpar.summary(pmod)
mean sd q0.025 q0.5 q0.975 mode
SD for the Gaussian observations 0.28651 0.055529 0.19794 0.27906 0.41494 0.26358
SD for site 0.46599 0.219836 0.18247 0.41706 1.02812 0.33867
bri.hyperpar.plot(pmod)
(lambda <- 3*sdres/tan(pi*0.99/2))
[1] 0.022066
halfcauchy <- "expression:
lambda = 0.022;
precision = exp(log_precision);
logdens = -1.5*log_precision-log(pi*lambda)-log(1+1/(precision*lambda^2));
log_jacobian = log_precision;
return(logdens+log_jacobian);"
hcprior <- list(prec = list(prior = halfcauchy))
formula <- nitrogen ~ f(site, model="iid", hyper = hcprior)
hmod <- inla(formula, family="gaussian", data=reeds)
bri.hyperpar.summary(hmod)
mean sd q0.025 q0.5 q0.975 mode
SD for the Gaussian observations 0.28763 0.056198 0.19809 0.28006 0.41769 0.26436
SD for site 0.42336 0.231601 0.14712 0.36566 1.03014 0.28080
reff <- pmod$marginals.random
x <- seq(-1.5,1.5,len=100)
d1 <- inla.dmarginal(x, reff$site[[1]])
d2 <- inla.dmarginal(x, reff$site[[2]])
d3 <- inla.dmarginal(x, reff$site[[3]])
rdf <- data.frame(nitrogen=x,density=c(d1,d2,d3),site=gl(3,100,labels=LETTERS[1:4]))
ggplot(rdf, aes(x=nitrogen, y=density, linetype=site))+geom_line()
bri.random.plot(pmod)
sdres <- sd(reeds$nitrogen)
pcprior <- list(prec = list(prior="pc.prec", param = c(3*sdres,0.01)))
formula <- nitrogen ~ f(site, model="iid", hyper = pcprior)
pmod <- inla(formula, family="gaussian", data=reeds, control.compute=list(config = TRUE))
psamp <- inla.posterior.sample(n=1000, pmod)
psamp[[1]]
$hyperpar
Precision for the Gaussian observations Precision for site
10.2559 6.2126
$latent
sample1
Predictor:01 2.90317
Predictor:02 2.89826
Predictor:03 2.90155
Predictor:04 2.90375
Predictor:05 2.90598
Predictor:06 2.65167
Predictor:07 2.64810
Predictor:08 2.65079
Predictor:09 2.64952
Predictor:10 2.65175
Predictor:11 3.46508
Predictor:12 3.46871
Predictor:13 3.47142
Predictor:14 3.46862
Predictor:15 3.47009
site:A 0.21266
site:B -0.03975
site:C 0.77777
(Intercept) 2.69090
$logdens
$logdens$hyperpar
[1] -5.4707
$logdens$latent
[1] 72.867
$logdens$joint
[1] 67.397
lvsamp <- t(sapply(psamp, function(x) x$latent))
colnames(lvsamp) <- row.names(psamp[[1]]$latent)
mean(lvsamp[,'site:C'] > lvsamp[,'site:A'])
[1] 0.992
Longitudinal Data
data(reading, package="brinla")
ggplot(reading, aes(agegrp, piat, group=id)) + geom_line()
library(lme4)
lmod <- lmer(piat ~ agegrp + (1|id), reading)
summary(lmod)
Linear mixed model fit by REML ['lmerMod']
Formula: piat ~ agegrp + (1 | id)
Data: reading
REML criterion at convergence: 1868.7
Scaled residuals:
Min 1Q Median 3Q Max
-1.946 -0.595 -0.028 0.507 2.868
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 29.8 5.46
Residual 44.9 6.70
Number of obs: 267, groups: id, 89
Fixed effects:
Estimate Std. Error t value
(Intercept) -11.538 2.249 -5.13
agegrp 5.031 0.251 20.04
Correlation of Fixed Effects:
(Intr)
agegrp -0.949
formula <- piat ~ agegrp + f(id, model="iid")
imod <- inla(formula, family="gaussian", data=reading)
imod$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -11.5352 2.26057 -15.9781 -11.5353 -7.0964 -11.5352 8.0785e-13
agegrp 5.0306 0.25287 4.5335 5.0306 5.5272 5.0306 1.0001e-12
bri.hyperpar.summary(imod)
mean sd q0.025 q0.5 q0.975 mode
SD for the Gaussian observations 6.7096 0.35619 6.0457 6.6959 7.4437 6.6650
SD for id 5.2988 0.62982 4.1389 5.2738 6.6074 5.2391
bri.hyperpar.plot(imod)
summary(imod$summary.random$id$mean)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-7.869 -2.934 -0.145 0.000 2.645 14.017
bri.random.plot(imod)
reading$cagegrp <- reading$agegrp - 8.5
lmod <- lmer(piat ~ cagegrp + (cagegrp|id), reading)
summary(lmod)
Linear mixed model fit by REML ['lmerMod']
Formula: piat ~ cagegrp + (cagegrp | id)
Data: reading
REML criterion at convergence: 1819.8
Scaled residuals:
Min 1Q Median 3Q Max
-2.654 -0.541 -0.150 0.385 3.281
Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 35.72 5.98
cagegrp 4.49 2.12 0.83
Residual 27.04 5.20
Number of obs: 267, groups: id, 89
Fixed effects:
Estimate Std. Error t value
(Intercept) 31.225 0.709 44.0
cagegrp 5.031 0.297 16.9
Correlation of Fixed Effects:
(Intr)
cagegrp 0.563
nid <- length(levels(reading$id))
reading$numid <- as.numeric(reading$id)
reading$slopeid <- reading$numid + nid
formula <- piat ~ cagegrp + f(numid, model="iid2d", n = 2*nid) + f(slopeid, cagegrp, copy="numid")
imod <- inla(formula, family="gaussian", data=reading)
imod$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 31.2241 0.69723 29.8520 31.2241 32.5947 31.2241 1.1108e-12
cagegrp 5.0305 0.29162 4.4567 5.0305 5.6036 5.0305 1.4099e-12
bri.hyperpar.summary(imod)
mean sd q0.025 q0.5 q0.975 mode
SD for the Gaussian observations 5.52166 0.305437 4.94884 5.51143 6.14773 5.49058
SD for numid (component 1) 5.68547 0.557750 4.66364 5.65933 6.85172 5.61227
SD for numid (component 2) 1.90716 0.265745 1.44226 1.88667 2.48399 1.84545
Rho1:2 for numid 0.94699 0.043666 0.82917 0.95948 0.98841 0.97515
postmean <- matrix(imod$summary.random$numid[,2],nid,2)
postmean <- sweep(postmean,2,imod$summary.fixed$mean,"+")
p <- ggplot(reading, aes(cagegrp, piat, group=id)) + geom_line(col=gray(0.95)) + xlab("centered age")
p+geom_abline(data=postmean,intercept=postmean[,1],slope=postmean[,2])
library(gridExtra)
sd.epsilon <- bri.hyper.sd(imod$internal.marginals.hyperpar[[1]],internal=TRUE)
sd.intercept <- bri.hyper.sd(imod$internal.marginals.hyperpar[[2]],internal=TRUE)
sd.slope <- bri.hyper.sd(imod$internal.marginals.hyperpar[[3]],internal=TRUE)
p1 <- ggplot(data.frame(sd.epsilon),aes(x,y))+geom_line()+ggtitle("Epsilon")+xlab("piat")+ylab("density")
p2 <- ggplot(data.frame(sd.intercept),aes(x,y))+geom_line()+ggtitle("Intercept")+xlab("piat")+ylab("density")
p3 <- ggplot(data.frame(sd.slope),aes(x,y))+geom_line()+ggtitle("Slope")+xlab("piat")+ylab("density")
p4 <- ggplot(data.frame(imod$marginals.hyperpar[[4]]),aes(x,y))+geom_line()+ggtitle("Rho")+ylab("density")
grid.arrange(p1,p2,p3,p4,ncol=2)
bri.hyperpar.plot(imod, together=FALSE)
formula <- log(piat) ~ cagegrp + f(numid, model="iid2d", n = 2*nid) + f(slopeid, cagegrp, copy="numid")
imod <- inla(formula, family="gaussian", data=reading)
bri.density.summary(imod$marginals.hyperpar[[4]])
mean sd q0.025 q0.5 q0.975 mode
0.059024 0.112174 -0.161646 0.059181 0.276275 0.059951
data(reading, package="brinla")
reading$id <- as.numeric(reading$id)
newsub <- data.frame(id=90, agegrp = c(6.5,8.5,10.5), piat=c(18, 25, NA))
nreading <- rbind(reading, newsub)
formula <- piat ~ agegrp + f(id, model="iid")
imod <- inla(formula, family="gaussian", data=nreading, control.predictor = list(compute=TRUE))
pm90 <- imod$marginals.fitted.values[[270]]
p1 <- ggplot(data.frame(pm90),aes(x,y))+geom_line()+xlim(c(20,60))
newsub=data.frame(id=90, agegrp = c(6.5,8.5,10.5), piat=c(NA, NA, NA))
nreading <- rbind(reading, newsub)
formula <- piat ~ agegrp + f(id, model="iid")
imodq <- inla(formula, family="gaussian", data=nreading, control.predictor = list(compute=TRUE))
qm90 <- imodq$marginals.fitted.values[[270]]
p1+geom_line(data=data.frame(qm90),aes(x,y),linetype=2)+xlab("PIAT")+ylab("density")
nsamp <- 10000
randprec <- inla.hyperpar.sample(nsamp, imod)[,1]
neweps <- rnorm(nsamp, mean=0, sd=1/sqrt(randprec))
newobs <- inla.rmarginal(nsamp, pm90) + neweps
dens <- density(newobs)
p1 + geom_line(data=data.frame(x=dens$x,y=dens$y),aes(x,y),linetype=2)
Classical Z-matrix Model
library(lme4)
mmod <- lmer(nitrogen ~ 1+(1|site), reeds)
Z <- getME(mmod, "Z")
X <- getME(mmod, "X")
library(INLA)
n <- nrow(reeds)
formula <- y ~ -1 + X + f(id.z, model="z", Z=Z)
imodZ <- inla(formula, data = list(y=reeds$nitrogen, id.z = 1:n, X=X))
library(brinla)
bri.hyperpar.summary(imodZ)
mean sd q0.025 q0.5 q0.975 mode
SD for the Gaussian observations 0.29079 0.05833 0.19817 0.28283 0.42608 0.26647
SD for id.z 0.31313 0.15740 0.10893 0.27909 0.71369 0.22175
imodZ$summary.random
$id.z
ID mean sd 0.025quant 0.5quant 0.975quant mode kld
1 1 -0.0041086 0.085306 -0.209637 0.00016443 0.156661 -0.00015857 0.0052617
2 2 -0.0041088 0.085306 -0.209637 0.00016426 0.156660 -0.00015874 0.0052617
3 3 -0.0041069 0.085306 -0.209633 0.00016583 0.156664 -0.00015718 0.0052616
4 4 -0.0041101 0.085306 -0.209639 0.00016323 0.156658 -0.00015975 0.0052617
5 5 -0.0041067 0.085306 -0.209633 0.00016596 0.156665 -0.00015705 0.0052616
6 6 -0.0504086 0.147214 -0.530901 -0.00406941 0.025792 -0.00097854 0.0010167
7 7 -0.0504092 0.147213 -0.530901 -0.00407005 0.025791 -0.00098011 0.0010176
8 8 -0.0504093 0.147213 -0.530901 -0.00407012 0.025791 -0.00098028 0.0010177
9 9 -0.0504095 0.147213 -0.530901 -0.00407031 0.025791 -0.00098072 0.0010179
10 10 -0.0504095 0.147213 -0.530901 -0.00407041 0.025791 -0.00098096 0.0010181
11 11 0.0533301 0.154185 -0.027216 0.00339512 0.557732 0.00069156 0.0015510
12 12 0.0533294 0.154185 -0.027216 0.00339500 0.557730 0.00069112 0.0015510
13 13 0.0533319 0.154186 -0.027216 0.00339549 0.557737 0.00069283 0.0015511
14 14 0.0533312 0.154186 -0.027216 0.00339536 0.557735 0.00069237 0.0015511
15 15 0.0533293 0.154184 -0.027216 0.00339496 0.557729 0.00069100 0.0015509
16 16 -0.0041052 0.085314 -0.209677 0.00018201 0.156701 -0.00016665 0.0052413
17 17 -0.0503191 0.147077 -0.530555 -0.00406629 0.025763 -0.00099760 0.0010821
18 18 0.0531454 0.153943 -0.027167 0.00337380 0.557112 0.00069253 0.0016927
data(meatspec, package="brinla")
trainmeat <- meatspec[1:172,]
testmeat <- meatspec[173:215,]
wavelengths <- seq(850, 1050, length=100)
modlm <- lm(fat ~ ., trainmeat)
rmse <- function(x,y) sqrt(mean((x-y)^2))
rmse(predict(modlm,testmeat), testmeat$fat)
[1] 3.814
plot(wavelengths,coef(modlm)[-1], type="l",ylab="LM Coefficients")
n <- nrow(meatspec)
X <- matrix(1,nrow = n, ncol= 1)
Z <- as.matrix(meatspec[,-101])
y <- meatspec$fat
y[173:215] <- NA
scaley <- 100
formula <- y ~ -1 + X + f(idx.Z, model="z", Z=Z)
zmod <- inla(formula, data = list(y=y/scaley, idx.Z = 1:n, X=X), control.predictor = list(compute=TRUE))
predb <- zmod$summary.fitted.values[173:215,1]*scaley
rmse(predb, testmeat$fat)
[1] 1.9027
rcoef <- zmod$summary.random$idx.Z[216:315,2]
plot(wavelengths, rcoef, type="l", ylab="Ridge Coefficients")
int.fixed <- list(prec = list(initial = log(1.0e-9), fixed=TRUE))
u.prec <- list(prec = list(param = c(1.0e-3, 1.0e-3)))
epsilon.prec <- list(prec = list(param = c(1.0e-3, 1.0e-3)))
idx.X <- c(1, rep(NA,100))
idx.Z <- c(NA, 1:100)
scaley <- 100
formula <- y ~ -1 + f(idx.X, model="iid", hyper = int.fixed) + f(idx.Z, model="iid", hyper = u.prec)
amod <- inla(formula, data = list(y=y/scaley, idx.X=idx.X, idx.Z=idx.Z), control.predictor = list(A=cbind(X, Z),compute=TRUE), control.family = list(hyper = epsilon.prec))
predb <- amod$summary.fitted.values[173:215,1]
rmse(predb, testmeat$fat/scaley)*scaley
[1] 1.9013
Generalized Linear Mixed Models
Poisson GLMM
data(nitrofen, package="boot")
head(nitrofen)
conc brood1 brood2 brood3 total
1 0 3 14 10 27
2 0 5 12 15 32
3 0 6 11 17 34
4 0 6 12 15 33
5 0 6 15 15 36
6 0 5 14 15 34
library(dplyr)
library(tidyr)
lnitrofen <- select(nitrofen, -total) %>% mutate(id=1:nrow(nitrofen)) %>% gather(brood,live,-conc,-id) %>% arrange(id)
lnitrofen$brood <- factor(lnitrofen$brood,labels=1:3)
head(lnitrofen)
conc id brood live
1 0 1 1 3
2 0 1 2 14
3 0 1 3 10
4 0 2 1 5
5 0 2 2 12
6 0 2 3 15
lnitrofen$jconc <- lnitrofen$conc + rep(c(-10,0,10),50)
ggplot(lnitrofen, aes(x=jconc,y=live, shape=brood)) + geom_point(position = position_jitter(w = 0, h = 0.5)) + xlab("Concentration")
library(lme4)
glmod <- glmer(live ~ I(conc/300)*brood + (1|id), nAGQ=25, family=poisson, data=lnitrofen)
summary(glmod)
Generalized linear mixed model fit by maximum likelihood (Adaptive Gauss-Hermite Quadrature, nAGQ = 25) ['glmerMod']
Family: poisson ( log )
Formula: live ~ I(conc/300) * brood + (1 | id)
Data: lnitrofen
AIC BIC logLik deviance df.resid
313.9 335.0 -150.0 299.9 143
Scaled residuals:
Min 1Q Median 3Q Max
-2.208 -0.606 -0.008 0.618 3.565
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.0911 0.302
Number of obs: 150, groups: id, 50
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.6386 0.1367 11.99 < 2e-16
I(conc/300) -0.0437 0.2193 -0.20 0.84
brood2 1.1688 0.1377 8.48 < 2e-16
brood3 1.3512 0.1351 10.00 < 2e-16
I(conc/300):brood2 -1.6730 0.2487 -6.73 1.7e-11
I(conc/300):brood3 -1.8312 0.2451 -7.47 7.9e-14
Correlation of Fixed Effects:
(Intr) I(c/300) brood2 brood3 I(/300):2
I(conc/300) -0.831
brood2 -0.700 0.576
brood3 -0.714 0.587 0.735
I(cn/300):2 0.529 -0.609 -0.804 -0.571
I(cn/300):3 0.538 -0.617 -0.569 -0.804 0.613
formula <- live ~ I(conc/300)*brood + f(id, model="iid")
imod <- inla(formula, family="poisson", data=lnitrofen)
imod$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 1.639578 0.13599 1.36778 1.641091 1.90290 1.644199 6.3039e-12
I(conc/300) -0.041354 0.21789 -0.47284 -0.040367 0.38417 -0.038459 2.2411e-12
brood2 1.163779 0.13764 0.89730 1.162479 1.43753 1.159878 1.2599e-13
brood3 1.345928 0.13505 1.08493 1.344500 1.61502 1.341638 1.7652e-13
I(conc/300):brood2 -1.663449 0.24847 -2.15557 -1.662005 -1.17976 -1.659093 1.7982e-15
I(conc/300):brood3 -1.820739 0.24499 -2.30621 -1.819229 -1.34406 -1.816188 9.2485e-15
bri.hyperpar.summary(imod)
mean sd q0.025 q0.5 q0.975 mode
SD for id 0.29335 0.057522 0.18809 0.29033 0.41501 0.28449
library(reshape2)
mf <- melt(imod$marginals.fixed)
cf <- spread(mf,Var2,value)
names(cf)[2] <- 'parameter'
ggplot(cf,aes(x=x,y=y))+geom_line()+facet_wrap(~ parameter, scales="free")+geom_vline(xintercept=0)+ylab("density")
bri.fixed.plot(imod)
bri.fixed.plot(imod,together=TRUE)
multeff <- exp(imod$summary.fixed$mean)
names(multeff) <- imod$names.fixed
multeff[-1]
I(conc/300) brood2 brood3 I(conc/300):brood2 I(conc/300):brood3
0.95949 3.20201 3.84175 0.18948 0.16191
sden <- data.frame(bri.hyper.sd(imod$marginals.hyperpar[[1]]))
ggplot(sden,aes(x,y)) + geom_line() + ylab("density") + xlab("linear predictor")
bri.hyperpar.plot(imod)
mden <- data.frame(inla.tmarginal(function(x) exp(1/sqrt(x)), imod$marginals.hyperpar[[1]]))
ggplot(mden,aes(x,y)) + geom_line() + ylab("density") + xlab("multiplicative")
formula <- live ~ I(conc/300)*brood + f(id, model="iid")
imod <- inla(formula, family="poisson", data=lnitrofen, control.compute=list(dic=TRUE))
formula <- live ~ log(conc+1)*brood + f(id, model="iid")
imod2 <- inla(formula, family="poisson", data=lnitrofen, control.compute=list(dic=TRUE))
c(imod$dic$dic, imod2$dic$dic)
[1] 786.11 842.06
mreff <- imod$summary.random$id$mean
qqnorm(mreff)
qqline(mreff)
formula <- live ~ I(conc/300)*brood + f(id, model="iid")
imod <- inla(formula, family="poisson", data=lnitrofen, control.compute=list(cpo=TRUE))
plot(log(imod$cpo$cpo),ylab="log(CPO)")
lnitrofen[which.min(imod$cpo$cpo),]
conc id brood live jconc
134 310 45 2 10 310
lnitrofen %>% filter(brood == 2, conc==310)
conc id brood live jconc
1 310 41 2 0 310
2 310 42 2 0 310
3 310 43 2 0 310
4 310 44 2 0 310
5 310 45 2 10 310
6 310 46 2 0 310
7 310 47 2 0 310
8 310 48 2 0 310
9 310 49 2 0 310
10 310 50 2 0 310
pit <- imod$cpo$pit
n <- length(pit)
uniquant <- (1:n)/(n+1)
logit <- function(p) log(p/(1-p))
plot(logit(uniquant), logit(sort(pit)), xlab="uniform quantiles", ylab="Sorted PIT values")
abline(0,1)
which.max(pit)
[1] 134
sdu <- 0.3
pcprior <- list(prec = list(prior="pc.prec", param = c(3*sdu,0.01)))
formula = live ~ I(conc/300)*brood + f(id, model="iid", hyper = pcprior)
imod2 = inla(formula, family="poisson", data=lnitrofen)
bri.hyperpar.summary(imod2)
mean sd q0.025 q0.5 q0.975 mode
SD for id 0.30966 0.056721 0.20718 0.30634 0.43046 0.30016
lnitrofen$obsid = 1:nrow(nitrofen)
formula <- live ~ I(conc/300)*brood + f(id, model="iid") + f(obsid, model="iid")
imodo <- inla(formula, family="poisson", data=lnitrofen)
bri.hyperpar.summary(imodo)
mean sd q0.025 q0.5 q0.975 mode
SD for id 0.293850 0.0564034 0.1951276 0.2899458 0.415677 0.283363
SD for obsid 0.010303 0.0062123 0.0037649 0.0085252 0.027042 0.006218
Binary GLMM
data(ohio, package="brinla")
table(ohio$smoke)/4
0 1
350 187
xtabs(resp ~ smoke + age, ohio)/c(350,187)
age
smoke -2 -1 0 1
0 0.16000 0.14857 0.14286 0.10571
1 0.16578 0.20856 0.18717 0.13904
library(lme4)
modagh <- glmer(resp ~ age + smoke + (1|id), nAGQ=25, family=binomial, data=ohio)
summary(modagh)
Generalized linear mixed model fit by maximum likelihood (Adaptive Gauss-Hermite Quadrature, nAGQ = 25) ['glmerMod']
Family: binomial ( logit )
Formula: resp ~ age + smoke + (1 | id)
Data: ohio
AIC BIC logLik deviance df.resid
1603.3 1626.0 -797.6 1595.3 2144
Scaled residuals:
Min 1Q Median 3Q Max
-1.373 -0.201 -0.177 -0.149 2.508
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 4.69 2.16
Number of obs: 2148, groups: id, 537
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.1015 0.2191 -14.16 <2e-16
age -0.1756 0.0677 -2.60 0.0095
smoke 0.3986 0.2731 1.46 0.1444
Correlation of Fixed Effects:
(Intr) age
age 0.244
smoke -0.493 -0.008
formula <- resp ~ age + smoke + f(id, model="iid")
imod <- inla(formula, family="binomial", data=ohio)
imod$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -2.99276 0.20434 -3.416832 -2.98475 -2.614898 -2.96902 2.1272e-14
age -0.16665 0.06283 -0.290606 -0.16646 -0.043901 -0.16606 1.8456e-14
smoke 0.39136 0.23954 -0.077759 0.39068 0.863758 0.38935 1.3287e-12
ilogit <- function(x) exp(x)/(1 + exp(x))
ilogit(imod$summary.fixed[1,c(3,4,5)])
0.025quant 0.5quant 0.975quant
(Intercept) 0.031774 0.04812 0.068186
exp(imod$summary.fixed[-1, c(3,4,5)])
0.025quant 0.5quant 0.975quant
age 0.74781 0.84666 0.95705
smoke 0.92519 1.47798 2.37206
bri.hyperpar.summary(imod)
mean sd q0.025 q0.5 q0.975 mode
SD for id 1.9258 0.16156 1.6253 1.9196 2.2603 1.9097
exp(bri.hyperpar.summary(imod)[3:5])
[1] 5.0800 6.8180 9.5858
table(xtabs(resp ~ id, ohio))
0 1 2 3 4
355 97 44 23 18
library(gridExtra)
p1 <- ggplot(data.frame(imod$marginals.fixed[[1]]),aes(x,y))+geom_line()+xlab("logit")+ylab("density")+ggtitle("Intercept")
p2 <- ggplot(data.frame(imod$marginals.fixed[[2]]),aes(x,y))+geom_line()+xlab("logit")+ylab("density")+ggtitle("age")
p3 <- ggplot(data.frame(imod$marginals.fixed[[3]]),aes(x,y))+geom_line()+xlab("logit")+ylab("density")+ggtitle("smoke")
sden <- data.frame(bri.hyper.sd(imod$marginals.hyperpar[[1]]))
p4 <- ggplot(sden,aes(x,y)) + geom_line() + xlab("logit") + ylab("density")+ggtitle("SD(u)")
grid.arrange(p1,p2,p3,p4,ncol=2)
inla.pmarginal(0,imod$marginals.fixed$smoke)
[1] 0.050916
formula <- resp ~ age + smoke + f(id, model="iid")
imod <- inla(formula, family="binomial", data=ohio, control.compute=list(dic=TRUE,waic=TRUE))
formula <- resp ~ factor(age) + smoke + f(id, model="iid")
imod1 <- inla(formula, family="binomial", data=ohio, control.compute=list(dic=TRUE,waic=TRUE))
formula <- resp ~ factor(age)*smoke + f(id, model="iid")
imod2 <- inla(formula, family="binomial", data=ohio, control.compute=list(dic=TRUE,waic=TRUE))
c(imod$dic$dic, imod1$dic$dic, imod2$dic$dic)
[1] 1466.3 1463.2 1463.5
c(imod$waic$waic, imod1$waic$waic, imod2$waic$waic)
[1] 1418.9 1416.3 1417.7
exp(imod1$summary.fixed[-1, c(3,4,5)])
0.025quant 0.5quant 0.975quant
factor(age)-1 0.74517 1.08375 1.57718
factor(age)0 0.65627 0.95965 1.40216
factor(age)1 0.38379 0.57802 0.86375
smoke 0.92488 1.48190 2.38575
ohio$obst = rep(1:4,537)
ohio$repl = ohio$id + 1
formula <- resp ~ factor(age) + smoke + f(obst, model="ar1", replicate = repl)
imod <- inla(formula, family="binomial", data=ohio)
exp(imod$summary.fixed[-1, c(3,4,5)])
0.025quant 0.5quant 0.975quant
factor(age)-1 0.74421 1.08461 1.58185
factor(age)0 0.65311 0.95868 1.40595
factor(age)1 0.37896 0.57416 0.86227
smoke 0.92547 1.48607 2.39799
bri.hyperpar.summary(imod)
mean sd q0.025 q0.5 q0.975 mode
SD for obst 1.97258 0.170482 1.66238 1.96343 2.33132 1.9443
Rho for obst 0.98471 0.014733 0.94575 0.98887 0.99856 0.9966
cmod <- inla(formula, family="binomial", data=ohio, control.inla = list(correct = TRUE))
cmod$summary.fixed
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -3.114881 0.28248 -3.708581 -3.101354 -2.59882 -3.075079 2.7827e-14
factor(age)-1 0.088258 0.19950 -0.302904 0.088034 0.48030 0.087609 1.5160e-13
factor(age)0 -0.045929 0.20463 -0.448273 -0.045802 0.35532 -0.045534 5.5515e-13
factor(age)1 -0.600099 0.22165 -1.041175 -0.598184 -0.17010 -0.594445 5.7201e-13
smoke 0.430558 0.26989 -0.096802 0.429313 0.96421 0.426899 1.7247e-12
bri.hyperpar.summary(cmod)
mean sd q0.025 q0.5 q0.975 mode
SD for obst 2.33088 0.218061 1.90502 2.33241 2.75712 2.34995
Rho for obst 0.96272 0.024311 0.90517 0.96735 0.99461 0.98362
sessionInfo()
R version 3.4.0 (2017-04-21)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.5
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] reshape2_1.4.2 tidyr_0.6.1 dplyr_0.5.0 gridExtra_2.2.1 ggplot2_2.2.1 brinla_0.1.0 INLA_17.06.30
[8] sp_1.2-4 lme4_1.1-13 Matrix_1.2-9 rmarkdown_1.4
loaded via a namespace (and not attached):
[1] Rcpp_0.12.10 compiler_3.4.0 nloptr_1.0.4 plyr_1.8.4 tools_3.4.0 digest_0.6.12
[7] evaluate_0.10 tibble_1.3.0 nlme_3.1-131 gtable_0.2.0 lattice_0.20-35 DBI_0.6-1
[13] yaml_2.1.14 stringr_1.2.0 knitr_1.15.1 rprojroot_1.2 grid_3.4.0 R6_2.2.0
[19] minqa_1.2.4 magrittr_1.5 backports_1.0.5 scales_0.4.1 htmltools_0.3.5 MASS_7.3-47
[25] splines_3.4.0 assertthat_0.2.0 colorspace_1.3-2 labeling_0.3 stringi_1.1.5 lazyeval_0.2.0
[31] munsell_0.4.3